---
title: "Workplace Mobility and Demographic by County"
author: ""
output: 
  flexdashboard::flex_dashboard:
    source_code: embed
---

```{r setup, include=FALSE}
library(highcharter)
library(flexdashboard)
library(ggplot2)
library(knitr)
library(kableExtra)


thm <- 
  hc_theme(
    colors = c("#1a6ecc", "#434348", "#90ed7d"),
    chart = list(
      backgroundColor = "transparent",
      style = list(fontFamily = "Source Sans Pro")
    ),
    xAxis = list(
      gridLineWidth = 1
    )
  )

```

Workplace Mobility
===================================== 

Column {data-width=500}
-----------------------------------------------------------------------
```{r}

# load ACS data
ACS_econ <- read.csv("../data/demographics/ACS_ECONOMIC_2018.csv",stringsAsFactors=FALSE)
ACS_social <- read.csv("../data/demographics/ACS_SOCIAL_2018.csv",stringsAsFactors=FALSE)

# rename GEO_ID to FIPS
names(ACS_econ)[1] <- "FIPS1"
names(ACS_social)[1] <- "FIPS2"

# extract abbreviated FIPS code
ACS_econ$FIPS1 = substr(ACS_econ$FIPS, 10, 14)
ACS_social$FIPS2 = substr(ACS_social$FIPS, 10, 14)

ACS_econ <- ACS_econ[-1,]
ACS_social <- ACS_social[-1,]


# laod SVI CDC data
SVI_dat <- read.csv("../data/demographics/SVI2018_US_COUNTY.csv")

# load mobility data
mobility <- read.csv("../data/mobility/county/county-data-wide-cleaned.csv")
mobility <- na.omit(mobility)

mobility$net_mob <- (mobility$X2020.03.22 + mobility$X2020.03.23 + mobility$X2020.03.24 + mobility$X2020.03.25 + mobility$X2020.03.26 + mobility$X2020.03.27 + mobility$X2020.03.28 + mobility$X2020.03.29)/8

# Aggregate all mobility for each county (average over categories to get net mobility)
dat = {}
mobility_agg <- data.frame(mobility$fips)
count = 1
  for (i in 6:ncol(mobility)){
    for (fip in mobility$fips){
    dat[count] <- sum(mobility[mobility$fips==fip,i])
    count = count + 1
    }
    dat <- dat/6
    mobility_agg <- cbind(mobility_agg,dat)
    names(mobility_agg)[i-4] <- names(mobility)[i]
    dat = {}
    count = 1
}

# clean up to get each county by fips
mobility_agg<-unique(mobility_agg)
names(mobility_agg)[1] <- "County_FIPS"


mobility_work <- cbind("FIPS3"=mobility$fips[mobility$Category=="Workplace"],"net_work_mob"=mobility$net_mob[mobility$Category=="Workplace"])

# select relevant columns from SVI data
SVI_sub <- SVI_dat[,c("FIPS", "STATE", "ST_ABBR","E_TOTPOP","EP_POV","EP_UNEMP","EP_PCI","EP_AGE65","EP_MINRTY","EP_MUNIT")]
ACS_social_sub <- ACS_social[,c("FIPS2","DP05_0037PE","DP05_0038PE","DP05_0039PE","DP05_0044PE","DP05_0052PE","DP05_0057PE","DP05_0071PE","DP05_0002PE","DP05_0003PE")]
ACS_econ_sub <- ACS_econ[,c("FIPS1","DP03_0002PE","DP03_0009PE","DP03_0025E","DP03_0027PE","DP03_0028PE","DP03_0029PE","DP03_0030E","DP03_0031PE","DP03_0062E","DP03_0063E")]


ACS_social_sub[,2:10] <- as.numeric(unlist(ACS_social_sub[,2:10]))
ACS_econ_sub <- na.omit(ACS_econ_sub)

ACS_econ_sub[,2:11] <- as.numeric(unlist(ACS_econ_sub[,2:11]))
ACS_econ_sub <- na.omit(ACS_econ_sub)


# join SVI, ACS, and mobility into one frame
full_dat <- merge(mobility_agg, SVI_sub, by.x="County_FIPS", by.y="FIPS")
full_dat <- merge(full_dat, ACS_social_sub,  by.x="County_FIPS", by.y="FIPS2")
full_dat <- merge(full_dat, ACS_econ_sub,  by.x="County_FIPS", by.y="FIPS1")
full_dat <- merge(full_dat, mobility_work,  by.x="County_FIPS", by.y="FIPS3")

# get average change in mobility for last week in data
full_dat$net_mob <- (full_dat$X2020.03.22 + full_dat$X2020.03.23 + full_dat$X2020.03.24 + full_dat$X2020.03.25 + full_dat$X2020.03.26 + full_dat$X2020.03.27 + full_dat$X2020.03.28 + full_dat$X2020.03.29)/8


```
### Workplace Mobility vs. Income

```{r}
ggplot(full_dat, aes(x=as.numeric(DP03_0062E),y=net_work_mob)) + geom_point(color="skyblue2") + geom_smooth(method='lm',color="skyblue4") + xlab("Median Household Income") + ylab("Average Change in Workplace Mobility 3/22-3/29")
```

### Workplace Mobility vs. Commute Time

```{r}
ggplot(full_dat, aes(x=as.numeric(DP03_0025E),y=net_work_mob)) + geom_point(color="aquamarine3") + geom_smooth(method='lm',color="aquamarine4") + xlab("Mean Commute Time (minutes)") + ylab("Average Change in Workplace Mobility 3/22-3/29")

```

Column {.tabset data-width=500}
-----------------------------------------------------------------------

### Sector Plots 

```{r}
options(repr.plot.width = 2, repr.plot.height = 1.5)
ggplot(full_dat, aes(x=as.numeric(DP03_0028PE),y=net_work_mob)) + geom_point(color="slateblue3") + geom_smooth(method='lm',color="slateblue4") + ggtitle("Service") + xlab("Percent Service Jobs") + ylab("Average Change in Workplace Mobility 3/22-3/29")

options(repr.plot.width = 2, repr.plot.height = 1.5)
ggplot(full_dat, aes(x=as.numeric(DP03_0029PE),y=net_work_mob)) + geom_point(color="slategray3") + geom_smooth(method='lm',color="slategray4") + ggtitle("Sales/Office") + xlab("Percent Sales/Office Jobs") + ylab("Average Change in Workplace Mobility 3/22-3/29")

options(repr.plot.width = 2, repr.plot.height = 1.5)
ggplot(full_dat, aes(x=as.numeric(DP03_0031PE),y=net_work_mob)) + geom_point(color="darkseagreen3") + geom_smooth(method='lm',color="darkseagreen4") + ggtitle("Production/Transport") + xlab("Percent Production/Transport Jobs") + ylab("Average Change in Workplace Mobility 3/22-3/29")

```

### Correlation Table

```{r}

# Income correlation
inc <- cor.test(full_dat$DP03_0062E,full_dat$net_work_mob)
# Commute correlation
com <- cor.test(full_dat$DP03_0025E,full_dat$net_work_mob)
# Service correlation
serv <-cor.test(full_dat$DP03_0028PE,full_dat$net_work_mob)
# Sales/Office 
sale <- cor.test(full_dat$DP03_0029PE,full_dat$net_work_mob)
# Production/Transport
prod <- cor.test(full_dat$DP03_0031PE,full_dat$net_work_mob)

correl <- c(-0.480,-0.176,0.024,-0.284,0.302)
p_values <- c("< 2.2e-16","< 2.2e-16","0.2605","< 2.2e-16","< 2.2e-16")

cor_tab <- cbind(correl, p_values)
rownames(cor_tab) <- c("Income","Commute","Service","Sales/Office","Production/Transport")
colnames(cor_tab) <- c("Pearson's r","p-value")

kable(cor_tab) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))

```

Examining the Pearson's r correlation coefficient for each variable and the average workplace mobility allows us to understand how the onset of COVID-19 has affected different aspects of citizens' work lives. As shown in the table above, both median income and commute time have significant negative correlations with workplace mobility. That is, the higher a person's income, the more they have decreased their mobility since the virus's spread. This is intuitive because higher-earning jobs are more likely to be successfully moved to a work-from-home format. This is further corroborated by the correlation breakdown by sector. Jobs that are sales and office-based have a significant negative correlation with mobility, indicating a large shift to remote work. Service jobs have shown no significant change in mobility, indicating they have remained working in the same fashion has before. Production and transportation jobs have shown a significant postive realtionship with mobility, indicating this sector has seen more of a boom with COVID. This basic demographic data gives valuable isight into how the workforce has collectively been affected by these unprecedent circumstances. 

Demographics
===================================== 

Column {data-width=500}
-----------------------------------------------------------------------
```{r}
ggplot(full_dat, aes(x=as.numeric(EP_POV),y=net_mob)) + geom_point(color="goldenrod2") + geom_smooth(method='lm',color="goldenrod4") + ggtitle("Persons Below Poverty Estimate") + xlab("Percent Below Poverty Line") + ylab("Average Change in Mobility 3/22-3/29")
```

```{r}
ggplot(full_dat, aes(x=as.numeric(DP05_0037PE),y=net_mob)) + geom_point(color="darkorange2") + geom_smooth(method='lm',color="darkorange4") + ggtitle("White Mobility") + xlab("Percent White") + ylab("Average Change in Mobility 3/22-3/29")
```

```{r}
ggplot(full_dat, aes(x=as.numeric(DP05_0044PE),y=net_mob)) + geom_point(color="tomato3") + geom_smooth(method='lm',color="tomato4") + ggtitle("Asian Mobility") + xlab("Percent Asian") + ylab("Average Change in Mobility 3/22-3/29")
```

Column {data-width=500}
-----------------------------------------------------------------------
```{r}
ggplot(full_dat, aes(x=as.numeric(EP_MINRTY),y=net_mob)) + geom_point(color="lightpink2") + geom_smooth(method='lm',color="lightpink4") + ggtitle("Minority Mobility") + xlab("Percent Minority") + ylab("Average Change in Mobility 3/22-3/29")
```

```{r}
ggplot(full_dat, aes(x=as.numeric(DP05_0038PE),y=net_mob)) + geom_point(color="lightsalmon2") + geom_smooth(method='lm',color="lightsalmon4") + ggtitle("African-American Mobility") + xlab("Percent African-American") + ylab("Average Change in Mobility 3/22-3/29")
```

```{r}
ggplot(full_dat, aes(x=as.numeric(DP05_0071PE),y=net_mob)) + geom_point(color="indianred2") + geom_smooth(method='lm',color="indianred4") + ggtitle("Hispanic/Latino Mobility") + xlab("Percent Hispanic/Latino") + ylab("Average Change in Mobility 3/22-3/29")
```

Mobility Trends by State
===================================== 

Column {data-width=500}
-----------------------------------------------------------------------

### Workplace Mobility vs. Time

```{r}
library(dplyr)
library(tidyr)
library(data.table)
library(plotly)
library(crosstalk)

# Aggregate all columns in c by weighted average, using E_TOTPOP as the weights
c <- colnames(full_dat)[!colnames(full_dat) %in% c("STATE", "E_TOTPOP", "County_FIPS")]
state_mobility <- data.table(full_dat, key = "STATE")
for(n in c){
 state_mobility[,
     (n) := weighted.mean(get(n), E_TOTPOP),
     #with = FALSE,
     by = STATE]
}
# Remove uninteresting columns
state_mobility <- state_mobility %>% select(-contains("DP")) %>% select(-contains("EP"))
state_mobility <- subset(state_mobility, select = -c(County_FIPS, E_TOTPOP, net_mob, net_work_mob, ST_ABBR ))
# Remove duplicate rows
state_mobility <- as.data.frame(state_mobility) %>% distinct()


state_policies <- read.csv("../data/policies/covid_us_state_policies.csv")
colnames(state_policies) <- toupper(colnames(state_policies))
state_policies$STATE <- toupper(state_policies$STATE)
state_policies <- subset(state_policies, select = c(STATE, STATE.OF.EMERGENCY,  STAY.AT.HOME..SHELTER.IN.PLACE))
state_policies$STATE.OF.EMERGENCY <- as.Date(state_policies$STATE.OF.EMERGENCY, format='%m/%d/%Y')
state_policies$STAY.AT.HOME..SHELTER.IN.PLACE <- as.Date(state_policies$STAY.AT.HOME..SHELTER.IN.PLACE, format='%m/%d/%Y')
state_policies[state_policies == 0] <- NA

# join the tables
state_mobility <- merge(x=state_mobility, y=state_policies, by="STATE", all.x=TRUE)

# Wide form to long form
state_mobility <- gather(state_mobility, date, mobility, X2020.02.17:X2020.03.29, factor_key=FALSE)
state_mobility$date <- as.Date(state_mobility$date, format='X%Y.%m.%d')

df.state<-group_by(state_mobility, STATE)
sd.state <- SharedData$new(df.state, ~STATE)
p0<-
   plot_ly(data=sd.state,x = ~date, y = ~mobility, height=750) %>% 
    add_lines(name="Selected",colors="gray",alpha=0.7) %>% 
    add_lines(name="State of Emergency Declared", x= ~STATE.OF.EMERGENCY) %>% 
  add_lines(name="Stay at Home Order Given", x= ~STAY.AT.HOME..SHELTER.IN.PLACE) %>% 
    add_lines(name="All states",data=state_mobility,x=~date,y=~mobility,
              colors="black",color=~STATE,alpha=0.1,showlegend=F,hoverinfo="none") %>%
     layout(title = "Mobility Trends by State",xaxis = list(title="Date"), 
            yaxis = list(title="% Change in Mobility from Baseline", range=c(-40, 20)))
bscols(widths=c(2,10),
  list(filter_select("state", "Select state to highlight for plot", sd.state, ~STATE,multiple=FALSE)),
  subplot(p0,titleY=T)
  )
```